%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

% ESTIMATING Full MODEL PARAMETERS			       		;
% T and F0 are now known at time=0                      ;
% T and F0 can now be correlated                        :
% Multi-Stage Optimization                              ;
% Oct 20, 2014 -- Paulina Oliva & Chris Severen         ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

clear all;
warning off all;
format compact
format shortg

load matdata_201312_samp1.mat

%*******************************************************;
% Execution and Optimization Options                    ;
%*******************************************************;

rp.parl    = 1;     % Run in parallel: 0 - serial, 1 - parallel
rp.cores   = 12;    % Only used if parallel==1
rp.k       = 1500;
rp.m       = 100;
rp.lam     = 0.5;     % Smoothing parameter
rp.tol     = 1e-15; % Tolerance in LL
rp.theseed = 25000; % RNG initilization
rp.cc      = 0.0001; % Convergence criterion
  
%*******************************************************;
% Extracting Data and Treatments                        ;
%*******************************************************;

obs     = size(data,1);

A       = data(:,1);
R       = data(:,2);
DP      = data(:,3);
DC_DP   = data(:,4);
TG      = data(:,5);
MON     = data(:,6);
I_N     = data(:,7:57);

%*******************************************************;
% Pre-defining useful inputs                            ;
%*******************************************************;

rp.df   = 0.6;   % Discount Factor
N       = 0:1:50;
Nbar    = 35;

if rp.parl==0
    % Note: If run serial, it is faster to generate everything once, as
    % done here. If run in parallel, it is faster to generate these locally
    % than using up transfer bandwidth. This creates everything and
    % collects it into the structure snPass.
    
    rng(rp.theseed);
    snPass.s       = rand(obs*(3*rp.k + rp.m),1);
    
    % These two lines create a (k x m x 51) matrix, where the first floor is a
    % (k x m) matrix of ones, the second a (k x m) matrix of twos, up to 51;
    % Nv is just N repeated k x 51 times;
    % the fourth and fifth are true false matrices for N vs. Nbar;
    Ntr            = permute(N,[3,1,2]);
    snPass.N3D     = Ntr(ones(rp.k,1),ones(rp.m,1),:);
    snPass.Nv      = N(ones(rp.k,1),:);
    snPass.Ngeq3D  = (snPass.N3D>=Nbar);
    snPass.N03D    = (snPass.N3D>0);
    snPass.Ngeqv   = (snPass.Nv>=Nbar);
    snPass.N0v     = (snPass.Nv>0);
elseif rp.parl==1
    snPass.s       = 0;
end

%*******************************************************;
% Initial values                                        ;
%*******************************************************;
initial     = zeros(1,8);

initial(1)  = 3.2;   %muT
initial(2)  = 1.2;   %sdT
initial(3)  = 0.7;    %rho

initial(4)  = 450;   %sdf0
initial(5)  = 100;   %sdf1

initial(6)  = 100;   %muF
initial(7)  = 0;   %ftreat
initial(8)  = -100;   %fmoni

%*******************************************************;
% Function call, start min                              ;
%*******************************************************;
    
if rp.parl==1
    matlabpool('open',rp.cores)
end

[theta_hat,value,iter_2s,diff,nfevaltot] = Full_3Stages(initial,DP,I_N,A,R,TG,MON,rp,snPass);
                   
if rp.parl==1
    matlabpool close
end
   
save results.mat
